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Abstract 

' In this report we review modern nonlinearity methods that can be used in the preterm 

, birth analysis. The nonhnear analysis of uterine contraction signals can provide information 

regarding physiological changes during the menstrual cycle and pregnancy. This information 
can be used both for the preterm birth prediction and the preterm labor control. 
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. 1 Introduction 

cr: 

Preterm birth is the single most important cause of perinatal mortahty in North America and 
Europe [Berkowitz Papiernik, 1993| . Developed countries represent 20 % of the population 
^ ' in the world, but only 12 % of human births annually, while 98 % of medical publications 

■ are issued from these areas. Reproductive patterns in the developed world, for the last three 

, decades, are different from elsewhere and during the first 70 years of the 20th century. A 

I major difference is in the number of children in families but also, and mainly, in the ages 

' at first pregnancies in primiparae (approaching now 30 years in many developed countries) 



O , [Robillard et al., 2007 . A British study [Rush et al., 1976] has estimated that preterm birth 



I accounts for 85 percent of early neonatal deaths that are not caused by lethal congenital mal- 

formations. In the United States, the smallest of the preterm infantsthose weighing less than 
. 750 g have been found to account for 41 percent of early neonatal deaths and 25 percent of 

r> I infant deaths [Overpeck et al., 1992 . Preterm birth is also a major determinant of neona- 



' tal and infant morbidity, including neuro-developmental handicaps, chronic respiratory prob- 

lems, intraventricular hemorrhage, infection, retrolental fibroplasia, and necrotizing enterocolitis 
Nat. Acad. Sci, 1985 Arias &: Tomich, 1993 . While there are inconsistent results with regard 



to long-term somatic growth among preterm infants Kitchen et al., 1992, Ross et al., 1990 , the 
risk of neurologic and developmental impairment during childhood is substantially elevated for 
the smallest survivors [Escobar et al., 1991 Veen et al, 1991] . In addition, the neonatal and 



long-term health care costs of preterm infants impose a considerable economic burden both on 
individual families and nationally [Nat. Acad. Sci, 1985| . 

While the frequency of births of low birth weight (less then 2,500 g) infants declined some- 
what in the United States between 1970 and 1980, this decline appears to have occurred primarily 
among full-term as opposed to preterm low birth weight infants [Kessel et al., 1984] . Further- 
more, although the infant mortality rate has decreased substantially since 1965, improvement in 
infant survival has been attributed principally to reduced birth weight-specific mortality rather 
than to changes in the birth weight or gestational age distribution [ Kleinman &: Kessel, 198T| . 
Despite the reduction in infant mortality, the rate in the United States remains considerably 
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higher than the rates in many other industriahzed countries. It is unhkely that there wih be 
further substantial improvement in infant survival in the United States unless a reduction in 
births of preterm low birth weight infants can be accomplished [Berkowitz fc: Papiernik, 1993 . 



Traditionally, prematurity was defined as a birth weight less than or equal to 2,500 g. How- 
ever, studies performed during the 1960s and 1970s revealed that this definition encompasses 
three distinct types of infants: those who are small because they were born too early, those 
who are small because their growth was retarded in utero, and those who are small because 
they were both premature and growth-retarded. The label prematurity has now been replaced 
by the terms low birth weight and preterm. According to current World Health Organization 
nomenclature [World Health Org. 1977] , low birth weight characterizes an infant who weighs 
less than 2,500 g (5 pounds and 8 ounces) at birth, and preterm refers to a birth that occurs at 
a gestational age of less than 37 completed weeks (less then 259 days). 

A cutoff of 37 weeks for defining preterm birth, albeit arbitrary, is now well established, but 
some studies have used cutoffs of 36 or 38 weeks of gestation. In addition, some investigators 
have only included preterm infants of a particular weight, such as those weighing less than 2,500 
g. Further subdivisions have also been used in the recent literature, such as very low birth 
weight to describe an infant weighing less than 1,500 g or 1,000 g at birth, and very preterm for 
gestational ages variously defined as less than 32, 33, or 34 weeks. Use of these lower cutoffs 
should be encouraged, since it will permit identification of risk factors that have the greatest 
impact on neonatal mortality [Berkowitz &: Papiernik, 1993| . 

For detecting contractions in uterine electromyography (EMG), fast numerical algorithms 
have been developed [Radharkrishnan et al., 2000[ . Recently, the analysis of uterine contraction 
signals has provided information regarding physiological changes during the menstrual cycle and 
pregnancy [Oczeretko et al, 2006[ . The authors have presented the cross- correlation and the 
wavelet cross- correlation methods to assess synchronization between contractions in different 
topographic regions of the uterus. It is important to identify time delays between uterine 
contractions, which may be of potential diagnostic significance in various pathologies. The cross- 
correlation was computed in a moving window with a width corresponding to approximately 
two or three contractions. As a result, the running cross-correlation function was obtained. 
The propagation % parameter assessed from this function allows quantitative description of 
synchronization in bivariate time series. In general, the uterine contraction signals are very 
complicated. Wavelet transforms provide insight into the structure of the time series at various 
frequencies (scales). To show the changes of the propagation % parameter along scales, a 
wavelet running cross-correlation was used. At first, the continuous wavelet transforms as the 
uterine contraction signals were received and afterwards, a running cross-correlation analysis 
was conducted for each pair of transformed time series. The findings of Oczeretko et al., 2006[ 
show that running functions are very useful in the analysis of uterine contractions. 

Fractal analysis serves another example for advanced data analysis to be utilized in the 
studies on uterine contractile activity. Fractals, a relatively new analytical concept of the 
last few decades, have been successfully applied in many areas of science and technology. 
Oczeretko et al., 2004[ performed a comparative study of ten methods of fractal analysis of 



signals of intrauterine pressure, and found significant differences between uterine contractions in 
healthy volunteers and women with primary dysmenorrhoea. There was a correlation between 
the adjacent elements of the investigated signals. Consequently, the values of fractal dimension 
can be objective measures for classifications of uterine contraction signals and may be used in 
studies on preterm labor. 

Nonlinear dynamics is one more example of advanced data analysis that could be imple- 
mented for uterine contractility signals [Pierzynski et al., 2007[ . In recent years the physiological 
signals obtained from the complex systems like the brain or the heart have been investigated for 
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possible deterministic chaotic behavior. The human uterus is undoubtedly a complex system. 
Smooth muscles comprising the myometrium interact in a complex manner. The techniques 
of surrogate data analysis have been used to testing for nonlinearity in the uterine contraction 
signals. The results showed that the spontaneous uterine contractions are considered to contain 
nonlinear features [Qczeretko et al., 2005| which indicated that nonlinear dynamics may increase 
the accessibility of data for the assessment of biology of uterus, and consequently, the nature of 
preterm labor. 



2 Techniques of Nonlinear Dynamics and Complex Data Anal- 
ysis 

In this section we review the concept of nonlinearity (as contrasted with the standard textbook 
linear algebra, analysis, statistics, etc.), together with its geometrical picture, and its extreme 
dynamics of chaotic behavior. 



2.1 Nonlinear Dynamics 

Recall that nonlinear dynamics is a modern language to talk about dynamical systems. It 
is a general theory of systems arising in physics, engineering, chemistry, biology, psychology, 
sociology and economics. Nonlinear dynamics has two extremes: at one end it is classical linear 
dynamics, fully predictable and controllable, as usually explained in engineering textbooks. On 
the other end, it is chaotic dynamics, or popularly, the "chaos theory. "0 Nonlinear dynamics 
includes both of these extremes, as well as all other natural and artificial dynamics. Its main 
keywords are (see, e.g., Ivancevic &: Ivancevic, 2006a| ): 



Dynamical system: A part of the world which can be seen as a self-contained entity 
with some temporal behavior. In nonlinear dynamics, speaking about a dynamical system 
usually means to speak about an abstract mathematical system which is a model for such 
an entity. Mathematically, a dynamical system is defined by its state and by its dynamics. 
A pendulum is an example for a dynamical system. 

State of a system: A number or a vector (i.e., a list of numbers) defining the state of 
the dynamical system uniquely. For the free (un-driven) pendulum, the state is uniquely 
defined by the angle 9 and the angular velocity = dO / dt. In the case of driving, the 
driving phase (j) is also needed because the pendulum becomes a non-autonomous system. 
In spatially extended systems, the state is often a field (a scalar-field or a vector-field). 
Mathematically spoken, fields are functions with space coordinates as independent vari- 
ables. The velocity vector-field of a fluid is a well-known example. 

Phase space: All possible states of the system. Each point in the phase-space corresponds 
to a unique state. In the case of the free pendulum, the phase-space has 2D whereas for 
driven pendulum it has 3D. The dimension of the phase-space is infinite in cases where 
the system state is defined by a field. 

Dynamics, or equation of motion: The causal relation between the present state and the 
next state in the future. It is a deterministic rule which tells us what happens in the next 
time step. In the case of a continuous time, the time step is infinitesimally small. Thus, 



^AU chaotic systems are nonlinear, but not all nonlinear systems are chaotic. Thus, chaotic dynamics is an 
extreme subset of nonlinear dynamics. 
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the equation of motion is an ordinary differential equation (ODE) (or a system of ODEs) : 

X = f{x), 

where X is the state and t is the time variable (overdot is the time derivative - as always). 
An example is the equation of motion of an un-driven and un-damped pendulum. In the 
case of a discrete time, the time steps are nonzero and the dynamics is a map: 

with the discrete time n. Note, that the corresponding physical time points i„ do not 
necessarily occur equidistantly. Only the order has to be the same. That is, 

n < m =^ tn < tm- 

The dynamics is linear if the causal relation between the present state and the next state 
is linear. Otherwise it is nonlinear. If we have the case in which the next state is not 
uniquely defined by the present one, this is generally an indication that the phase-space 
is not complete. Thus, there are important variables determining the state which had 
been forgotten. This is a crucial point while modelling a real-life systems. Beside this, 
there arc two important classes of systems where the phase-space is incomplete: the 
non-autonomuous and stochastic systems. A non-autonomous system has an equation of 
motion which depends explicitly on time. Thus, the dynamical rule governing the next 
state not only depends on the present state but also at the time it applies. A driven 
pendulum is a classical example of a non-autonomuous system. Fortunately, there is an 
easy way to make the phase-space complete: we simply include the time into the definition 
of the state. Mathematically, this is done by introducing a new state variable: t. Its 
dynamics reads 

t = 1, or tn+i = tji, 

depending on whether time is continuous or discrete. For the periodically driven pendula, 
it is also natural to take the driving phase as the new state variable. Its equation of motion 
reads 

9 = 2Trw, 

where w is the driving frequency (so that the angular driving frequency is 2Trw). On the 
other hand, in a stochastic system, the number and the nature of the variables necessary 

to complete the phase-space is usually unknown. Therefore, the next state can not be 
deduced from the present one. The deterministic rule is replaced by a stochastic one. 
Instead of the next state, it gives only the probabilities of all points in the phase-space to 
be the next state. 

• Orbit or trajectory: A solution of the equation of motion. In the case of continuous time, 
it is a curve in phase-space parameterized by the time variable. For a discrete system it 
is an ordered set of points in the phase-space. 

• Phase Flow: The mapping (or, map) of the whole phase-space of a continuous dynamical 
system onto itself for a given time step t. If t is an infinitesimal time step dt, the fiow is 
just given by the right-hand side of the equation of motion (i.e., /). In general, the fiow 
for a finite time step is not known analytically because this would be equivalent to have a 
solution of the equation of motion. 
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More precisely, a mathematical term dynamical system geometrically represents a vector- 
field in the system's phase-space manifold M (see |Ivancevic &: Ivancevic, 2006b| ), which upon 
integration (governed by the celebrated existence & uniqueness theorems for ordinary Differential 
Equations) defines a phase-flow in M (see Figured]). This phase-flow ft G M, describing the 
complete behavior of a dynamical system at every time instant, can be either linear, nonlinear 
or chaotic. 




Figure 1: Action of the phase-flow ft in the phase-space manifold M: (a) Trajec- 
tory of a single initial point x{t) G M , (b) Transporting the whole manifold M (see 
Ivancevic &: Ivancevic, 2006a| ). 

Before the advent of fast computers, solving a dynamical system required sophisticated 
mathematical techniques and could only be accomplished for a small class of linear dynamical 
systems. Numerical methods executed on computers have simplified the task of determining the 
orbits of a dynamical system. 

For simple dynamical systems, knowing the trajectory is often sufficient, but most dynamical 
systems are too complicated to be understood in terms of individual trajectories. The difficulties 
arise because: 

1. The systems studied may only be known approximately-the parameters of the system may 
not be known precisely or terms may be missing from the equations. The approximations used 
bring into question the validity or relevance of numerical solutions. To address these questions 
several notions of stability have been introduced in the study of dynamical systems, such as 
Lyapunov stability or structural stability. The stability of the dynamical system implies that 
there is a class of models or initial conditions for which the trajectories would be equivalent. 
The operation for comparing orbits to establish their equivalence changes with the different 
notions of stability. 

2. The type of trajectory may be more important than one particular trajectory. Some 
trajectories may be periodic, whereas others may wander through many different states of the 
system. Applications often require enumerating these classes or maintaining the system within 
one class. Classifying all possible trajectories has led to the qualitative study of dynamical 
systems, that is, properties that do not change under coordinate changes. Linear dynamical 
systems and systems that have two numbers describing a state are examples of dynamical systems 
where the possible classes of orbits are understood. 

3. The behavior of trajectories as a function of a parameter may be what is needed for an 
application. As a parameter is varied, the dynamical systems may have bifurcation points where 
the qualitative behavior of the dynamical system changes. For example, it may go from having 
only periodic motions to apparently erratic behavior, as in the transition to turbulence of a fluid. 

4. The trajectories of the system may appear erratic, as if random. In these cases it may 
be necessary to compute averages using one very long trajectory or many different trajectories. 
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The averages are well defined for ergodic systems and a more detailed understanding has been 
worked out for hyperbolic systems. Understanding the probabilistic aspects of dynamical systems 
has helped establish the foundations of statistical mechanics and of chaos. 



2.2 Geometrical View on Nonlinear Dynamics 

Geometrically speaking (see [Ivancevic &: Ivancevic, 2006b] ), linear systems (i.e., systems de- 
fined by linear differential equations) live in Euclidean spaces M'^, where is the dimension 
of the system. For their analysis the tools of linear algebra and calculus are used. The basic 
measure here is the global Euclidean metric form (or, Euclidean distance function) function on 
the system coordinates, defined by 



S 



\ 



N N 

^ {xi - yif, or 5^ = ^ {xi - yif 

i=l 1=1 



The more realistic nonlinear systems (i.e., systems defined by non-linear differential equa- 
tions) can be considered as local deformations of the closest linear systems and they live in 
smooth manifolds . For their analysis the tools of Riemannian geometry are used, with the 
local Riemannian metric form., defined by 

N N 

ds"^ = ''^^''^^ gijdxidxj , or ds^ = gijdx^dx^, (1) 
i=i j=i 

where gij is the Riemannian metric tensor, dx^ are differentials of the local coordinates, and 
the Einstein's summation convention (summing upon repeated indices if one is superscript- 
contravariant and the other is subscript-covariant) is in place. Besides giving the local distances 
between the points on the smooth manifold M^, the Riemannian metric form ([T|) defines the 
system's kinetic energy 

KE = —gijx''x\ (overdot means time derivative), 

giving the Euler-Lagrangian equations of motion in a geometrical form (for s = t) of geodesic 
equations 

d^x* dx^ dx^ , , 

where F*;!, are the so-called Christoffel symbols of the affine Levi-Civita connection of the man- 
ifold (see, e.g., [Ivancevic &: Ivancevic, 2007b| ). 

In the geometrical framework, the (in)stability of the trajectories is the (in)stability of the 
geodesies, and it is completely determined by the curvature properties of the underlying manifold 
according to the Jacobi equation of geodesic deviation 



Ivancevic & Ivancevic, 2007b 



D^J' i dx^ .dx^ 



+ R\ikm-r-J —r- - 0; (3) 



whose solution J, usually called Jacobi variation field, locally measures the distance between 
nearby geodesies; D/ds stands for the covariant derivative along a geodesic and R^ji^^y^ are the 
components of the Riemann curvature tensor. 
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In the special case of the Eisenhart metric [Eisenhart, 1929 | on an enlarged configuration 
space-time {{x^ = t,x^, . . . ,x^}, the Jacobi equation ([3]) reduces to the tangent dynamics 
equation [Ivancevic Sz Ivancevic, 2007b| 

2.3 Extreme Nonlinearity: Chaotic Dynamics 

On the other hand, chaotic dynamics, the most extreme case of nonlinear dynamics, is highly 
sensitive to both initial conditions (popularly referred to as the butterfly effect) and system 
parameters. As a result of this sensitivity, which manifests itself as an exponential growth of 
perturbations in the initial conditions, the behavior of chaotic systems appears to be random. 
This happens even though these systems are deterministic, meaning that their future dynamics 
are fully defined by their initial conditions, with no random elements involved. This behavior 
is known as deterministic chaos (see, e.g., [Ivancevic fc Ivancevic, 2006a| ). Chaotic behavior 
has been observed in the laboratory in a variety of systems including electrical circuits, lasers, 
oscillating chemical reactions, fluid dynamics, and mechanical and magneto-mechanical devices. 

Chaotic behavior is surprisingly complex and may arise in a simple noise-free system with 
few degrees of freedom. Chaos is globally stable and locally unstable. A chaotic system neither 
reaches a certain steady state nor repeats itself periodically. Thus, it is a pure disordered process 
since no points or patterns of points ever recur. However, its activity remains constrained within 
a certain boundary, which may suggest a new kind of order. Chaotic behavior resembles random 
noise, but is predictable in the short-term. This short-term predictability is useful in various 
domains ranging from weather forecasting to economic forecasting [Hilborn, 1994 . 



The most common route to chaos is a sequence of period-doubling bifurcations, caused 
by particular values of a nonlinear system parameter. The cases of most interest arise when 
the chaotic behavior takes place on an attractor, since then a large set of initial conditions 
will lead to orbits that converge to this chaotic region. In particular, the so-called strange 
attractors typically have a fractal structure. In all chaotic systems all particle trajectories 
diverge exponentially from one another, with a positive Lyapunov exponent. Besides, chaotic 
behavior is related to the non- equilibrium phase transition^ caused by system's topology chang^ 



Caiani et al., 1997 Franzosi et al., 1999, Franzosi et at, 2000J. 



■^Phase transitions (PT) are phenomena wliich bring about qualitative physical changes at the macro- 
scopic level in presence of the same microscopic forces acting among the constituents of a system (e.g., 
Solid ^ Liquid ^ Gas ^ Plasma). Their mathematical description requires to translate into quantitative 
terms the mentioned qualitative changes. The standard way of doing this is to consider how the values of thermo- 
dynamic observables, obtained in laboratory experiments, vary with temperature, or volume, or an external field, 
and then to associate the experimentally observed discontinuities at a PT to the appearance of some kind of sin- 
gularity entailing a loss of analyticity [Franzosi fc Pettini, 2004] , The so-called non-equilibrium phase transitions 
have been elaborated in the realm of synergeties 'Haken, 1983 , Hak en, 1993[ [Ivancevic &: Ivancev ic, 2006c as a 
ubiquitous route to self-organization in complex systems of various nature (e.g., in brain's function |Haken, 2002 ). 

^ Topology change means loss of the system's diffeomorphicity. Namely, the so-called topological theorem 
[Franzosi fc Pettini, 2004] says that non-analyticity is the "shadow" of a more fundamental phenomenon occurring 
in the system's configuration space: a topology change within the family of equipotential hyper- surfaces 

E„ = {(a;i, . . . ,xn) £ V{xx, . . . ,xn) = v}, 

where V — V{x) is the microscopic interaction potential expressed in the coordinates Xi. The largest Lyapunov 
exponent Ai (see subsection 12. 3.21 below) is computed by solving the above tangent dynamics equation ((J) so that 
we get [Franzosi fc Pettini, 2004[ 

Al = lim l/2tlog(EiIi[ji^(t) + Jf{t)]/Etilj?{0) + J?m). 
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Techniques have emerged to harness chaos to the benefit of humans. The chaotic behavior of 
a system may be artificiahy weakened or suppressed if it is undesirable. This concept is known as 
control of chaos. The first and the most important method of chaos control is the so-called OGY- 
method, developed by |Ott et al, 1990| . However, in recent years, a non-traditional concept of 
anti-control of chaos has emerged. Here, the non-chaotic dynamical system is transformed into 
a chaotic one by small controlled perturbation so that useful properties of a chaotic system can 
be exploited. This non-traditional concept has found applications in time and energy critical 
control applications such as navigation in a multi-body planetary system, fluid mixing and to 
secure information processing Chen &: Dong, 1998| . 



2.3.1 Lorenz Attractor 

Recall that a dynamical system may be defined as a deterministic rule for the time evo- 
lution of state observables. Well known examples are ODEs in which time is continuous 
[Ivancevic &: Ivancevic, 2006a| 



x(t) = f(x(t)), (x,f GM"); 
and iterative maps in which time is discrete: 

x(t + l) =g(x(t)), (x,gGE"). 



(5) 



(6) 



In the case of maps, the evolution law is straightforward: from x(0) one computes x(l), and then 
x(2) and so on. For ODE's, under rather general assumptions on f , from an initial condition 
x(0) one has a unique trajectory x(t) for t > Qtt, 1993] . Examples of regular behaviors (e.g., 
stable fixed-points, limit cycles) are well known, see Figure [2l 




-1 ^I.S 



Figure 2: Examples of regular attractors: fixed-point (left) and limit cycle (right). Note that 
limit cycles exist only in nonlinear dynamics. 

A rather natural question is the possible existence of less regular behaviors i.e., different 
from stable fixed-points, periodic or quasi-periodic motion. 

After the seminal works of Poincare, Lorenz, Smale, May, and Henon (to cite only the most 
eminent ones) it is now well established that the so called chaotic behavior is ubiquitous. As a 
relevant system, originated in the geophysical context, we mention the celebrated Lorenz system 
Lorenz, 1963[ [Sparrow, 1982 



X = —a{x — y) 
y = —xz + rx — y 
z = xy — bz 



(7) 



This system is related to the Rayleigh-Benard convection under very crude approximations. 
The quantity x is proportional the circulatory fiuid particle velocity; the quantities y and z are 
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related to the temperature profile; cr, b and r are dimensionless parameters. Lorenz studied the 
case with cr = 10 and 6 = 8/3 at varying r (which is proportional to the Rayleigh number). It 
is easy to see by linear analysis that the fixed-point (0,0,0) is stable for r < 1. For r > 1 it 
becomes unstable and two new fixed-points appear 



i±^/b{r-l), ±^/b{r-l), r - 1) 



(8) 



these are stable for r < Vc = 24.74. A nontrivial behavior, i.e., non periodic, is present for 
r > Tc, as is shown in Figure El 

In this 'strange', chaotic regime one has the so called sensitive dependence on initial con- 
ditions. Consider two trajectories, x(t) and x'(t), initially very close and denote with A(t) = 
||x'(t) — x(t)|| their separation. Chaotic behavior means that if A(0) — > 0, then as t — > oo one 
has A(t) 



A(0)expAit, with Ai > [Boffetta et al, 2001| . 
Let us notice that, because of its chaotic behavior and its dissipative nature, i.e., 

dx dii dz 
ox ay oz 



(9) 



the attractor of the Lorenz system cannot be a smooth surface. Indeed the attractor has a 
self-similar structure with a fractal dimension between 2 and 3. The Lorenz model (which 
had an important historical relevance in the development of chaos theory) is now considered a 
paradigmatic example of a chaotic system. 




Figure 3: Example of an aperiodic signal: the x variable of the Lorenz system d?]) as function 
of time t, for r = 28. 



2.3.2 Lyapunov Exponents 

The sensitive dependence on the initial conditions can be formalized in order to give it a quanti- 
tative characterization. The main growth rate of trajectory separation is measured by the first 
(or maximum) Lyapunov exponent, defined as (see, e.g., [Boffetta et al, 2001| ) 

Al = lim lim -In^^, (10) 
t^ooA(0)^Ot A(0) ^ ^ 

As long as A{t) remains sufficiently small (i.e., infinitesimal, strictly speaking), one can regard 
the separation as a tangent vector z{t) whose time evolution is 



dxj 



Ut) ■ Zj, (11) 
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and, therefore, 

Ai=limiln^M. (12) 

t^oot ||z(0)|| ^ ^ 

In principle, Ai may depend on the initial condition x(0), but this dependence disappears for 
ergodic systems. In general there exist as many Lyapunov exponents, conventionally written 
in decreasing order Ai > A2 > A3 > as the independent coordinates of the phase-space 
[Benettin et al., 1980 . Without entering the details, one can define the sum of the first k 



Lyapunov exponents as the growth rate of an infinitesimal /cD volume in the phase-space. In 
particular, Ai is the growth rate of material lines, Ai + A2 is the growth rate of 2D surfaces, and 



so on. A numerical widely used efficient method is due to Benettin et al., 1980 . 

It must be observed that, after a transient, the growth rate of any generic small perturba- 
tion (i.e., distance between two initially close trajectories) is measured by the first (maximum) 
Lyapunov exponent Ai, and Ai > means chaos. In such a case, the state of the system is 
unpredictable on long times. Indeed, if we want to predict the state with a certain tolerance A 
then our forecast cannot be pushed over a certain time interval Tp, called predictability time, 
given by IBoffetta et al, 2001| : 

The above relation shows that Tp is basically determined by 1/Ai, seen its weak dependence on 
the ratio A/A(0). To be precise one must state that, for a series of reasons, relation is too 



simple to be of actual relevance Boffetta et al., 2002 



2.3.3 Kolmogorov Sinai Entropy 

Deterministic chaotic systems, because of their irregular behavior, have many aspects in common 
with stochastic processes. The idea of using stochastic processes to mimic chaotic behavior, 
therefore, is rather natural [Chirikov, 1979 Benettin, 1984 . One of the most relevant and 



successful approaches is symbolic dynamics [Beck &: Schlogl, 1993] . For the sake of simplicity 
let us consider a discrete time dynamical system. One can introduce a partition A of the phase- 
space formed by N disjoint sets Ai, ...,A]sf. From any initial condition one has a trajectory 

x(0) ^x(l),x(2),...,x(n),... (14) 



dependently on the partition element visited, the trajectory (|T4j) . is associated to a symbolic 
sequence 

x(0) ^ -^1,^2, •••,^n, ••• (15) 

where in {n = 1,2,...,N) means that x(n) G Ai^ at the step n, for n = 1,2,.... The coarse- 
grained properties of chaotic trajectories are therefore studied through the discrete time process 

m- 

An important characterization of symbolic dynamics is given by the Kolmogorov-Sinai en- 
tropy (KS), defined as follows. Let C„ = (ii,i2, ■■■,in) be a generic 'word' of size n and P{Cn) 
its occurrence probability, the quantity [Boffetta et al., 2001 



Hn = snp[-^P{Cn)lnP{Cn)], (16) 



is called block entropy of the n— sequences, and it is computed by taking the largest value over all 
possible partitions. In the limit of infinitely long sequences, the asymptotic entropy increment 

hKS = lim Hn+i - Hn, (17) 
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is the Kolmogorov-Sinai entropy. The difference —Hn has the intuitive meaning of average 
information gain supphed by the (n + 1)— th symbol, provided that the previous n symbols are 
known. KS-entropy has an important connection with the positive Lyapunov exponents of the 
system |Utt, 1993| : 

hKS = Yl (18) 

Ai>0 

In particular, for low-dimensional chaotic systems for which only one Lyapunov exponent is 
positive, one has hxs = Ai. 

We observe that in ([16]) there is a technical difficulty, i.e., taking the sup over all the possible 
partitions. However, sometimes there exits a special partition, called generating partition, for 
which one finds that coincides with its superior bound. Unfortunately the generating parti- 
tion is often hard to find, even admitting that it exist. Nevertheless, given a certain partition, 
chosen by physical intuition, the statistical properties of the related symbol sequences can give 
information on the dynamical system beneath. For example, if the probability of observing a 
symbol (state) depends only by the knowledge of the immediately preceding symbol, the sym- 
bolic process becomes a Markov chain (see [Ivancevic fc: Ivancevic, 2006b| ) and all the statistical 
properties are determined by the transition matrix elements Wij giving the probability of ob- 
serving a transition i ^ j in one time step. If the memory of the system extends far beyond 
the time step between two consecutive symbols, and the occurrence probability of a symbol 
depends on k preceding steps, the process is called Markov process of order k and, in principle, 
a k rank tensor would be required to describe the dynamical system with good accuracy. It is 
possible to demonstrate that if Hn+i — Hn = hxs for n > k + 1, k is the (minimum) order of 
the required Markov process. It has to be pointed out, however, that to know the order of the 
suitable Markov process we need is of no practical utility if S> 1. 



2.4 Nonlinear Analysis of Time Series Data 

In this subsection, following Sharma, 2006| , we apply two basic techniques of nonlinear dy- 
namics to the heart inter-beat interval (IBI) time series data. The volunteer participant was a 
twenty year-old male pilot with 400 hours total flight experience and a commercial pilot license. 
He flew a single-engine aircraft in a fixed-base generic fiight simulator. A Polar Heart-rate Re- 
ceiver/Transmitter was used to record his heart IBI. The task of the participant in the flight 
segments Ml and M2 was to maintain manually a set airspeed, altitude, and heading as closely as 
possible. The only difference between the two segments. Ml and M2, was the level of perceived 
risk of mid-air collision involved if the participant was unable to execute the task adequately. 
First, a sample participant's heart IBI data at segments Ml and M2 are represented graphically 
and analyzed in terms of phase plots. Next, the measure of approximate entropy is employed to 
quantify the dynamics for each time series [Sharma, 2006| . 



2.4.1 Phase Plots 

The qualitative dynamics of a system may be represented in a graphical form. The Phase plot is 
the simplest graphical representation of the temporal information contained in a time series. It 
relates the current value of the time series to its preceding value. However, the preceding value 
is expressed in such a manner that it is approximately equal to the derivative of the current 
value. The phase plot provides a spatial representation of the evolving dynamics of a nonlinear 
system, giving some information on how the process evolves over time [Williams, 1997[ . 

Figure m represents two-dimensional phase plots employed to depict the graphical dynamics 
of a sample participant's heart IBI at fiight segments Ml and M2. The figure shows that the two 
plots are qualitatively different from each other in terms of their dynamics within the manifold. 
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Figure 4: Two-dimensional phase plots of heart inter-beat interval (IBI) for a pilot at the 
simulated flying segments Ml (left) and M2 (right; adapted from Sharma, 2006| ). 



which is represented by their trajectories. This suggests that they have different vector fields that 
govern their evolution within their manifolds. Further, a thorough analysis indicates that the 
system entropy (information content) at segment M2 is slightly higher than that at segment Ml. 
In other words, the heart IBI dynamics at segment M2 appears slightly less structured compared 
to that at segment Ml. However, some quantitative evidence may be required to support this 
observation, which has been done in the next section. The above observations suggest that heart 



IBI dynamics depends on the level of the pilot's perceived risk Sharma, 2006 



2.4.2 Approximate Entropy 

In general, entropy is concerned with information about the state of a dynamical system. Ap- 
proximate entropy (ApEn) is a useful general measure of nonlinear dynamics. The technical 
details of the ApEn index have been provided in Pincus, 1995| . ApEn is a measure quantifying 



the regularity of the time series [Pikkujamsa et al., 19991 . ApE n is defined as 

ApEn (m, r) = lim [0 (m, r, A^) — cp (m + 1, r, N)] 

N^oo 

where m is the maximum run length, r is the tolerance, is the total number of observa- 
tions in the time series, and (p (m, r, N) is the average value of the logarithm of proportion 
of vectors that are closer than the rover all vectors. ApEn has a relatively lower value for a 
structured/deterministic data set and a higher value for a disordered/less-deterministic data set 
|Ho et al., 1997] . 

An algorithm has been provided in [Kaplan et al., 1991 to compute ApEn. The algorithm 



detects similar patterns in a time series and then estimates the logarithmic likelihood that 
the following observation will differ from the previously observed pattern. Further, the pat- 
tern length and the pattern similarity criterion can be manipulated while running the program 
[Ho et al., 1997| . ApEn can be computed using no more than 100 observations [Pincus, 1995| . 

[Sharma, 2006[ employed ApEn to analyze the heart IBI dynamics of a pilot during flight 
segments Ml and M2. During the analysis, the value for the maximum run length (m) was 
set equal to 2 and the tolerance (r) was selected as 0.15 of the standard deviation as suggested 
in [Pincus, 1995[ . The resulting ApEn values at segments Ml and M2 were 0.452 and 0.551 
respectively. These quantitative results indicate that the heart IBI ApEn becomes higher at 
segment M2 compared to that at segment Ml. These results re-confirm the visual analysis of 
the heart IBI Phase Plot data in the previous section, where the entropy at segment M2 was 
believed to be slightly higher than that at segment Ml. This result may further imply that the 
structure of the heart IBI data becomes more disordered and unstructured during a higher level 
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of perceived-risk. Finally, these results suggest that the heart IBI ApEn may be a sensitive 
index of the level of perceived risk. 



2.5 Lyapunov Exponents in Time Series Data 

Recall that chaos arises from the exponential growth of infinitesimal perturbations, together 
with global folding mechanisms to guarantee boundedness of the solutions. This exponential 
instability is characterized by the spectrum of Lyapunov exponents 

[Eckmann &: Ruelle, 1985] . If one assumes a local decomposition of the phase space into direc- 
tions with different stretching or contraction rates, then the spectrum of exponents is the proper 
average of these local rates over the whole invariant set, and thus consists of as many exponents 
as there are space directions. The most prominent problem in time series analysis is that the 
physical phase space is unknown, and that instead the spectrum is computed in some embedding 
space. Thus the number of exponents depends on the reconstruction, and might be larger than 
in the physical phase space. Such additional exponents are called spurious, and there are several 
suggestions to either avoid them, or to identify them. Moreover, it is plausible that only as many 
exponents can be determined from a time series as are entering the Kaplan Yorke formula (see 
below). To give a simple example: Consider motion of a high-dimensional system on a stable 
limit cycle. The data cannot contain any information about the stability of this orbit against 
perturbations, as long as they are exactly on the limit cycle. For transients, the situation can 
be different, but then data are not distributed according to an invariant measure and the nu- 
merical values are thus difficult to interpret. Apart from these difficulties, there is one relevant 
positive feature: Lyapunov exponents are invariant under smooth transformations and are thus 
independent of the measurement function or the embedding procedure. They carry a dimension 
of an inverse time and have to be normalized to the sampling interval [Hegger et al, 1999 . 



2.5.1 The Maximal Lyapunov Exponent 

The maximal Lyapunov exponent can be determined without the explicit construction of a 
model for the time series. A reliable characterization requires that the independence of embed- 
ding parameters and the exponential law for the growth of distances are checked [Kantz, 1994 



Rosenstein et. al, 1993| explicitly. Consider the representation of the time series data as a trajec- 



tory in the embedding space, and assume that you observe a very close return s^' to a previously 



visited point Sn- Then one can consider the distance [Hegger et al, 1999 

Ao = Sn — Sn' 

as a small perturbation, which should grow exponentially in time. Its future can be read from 
the time series: 

A; = Sn+l — Sn'+l- 

If one finds that |A;| ~ Aqb^' then A is (with probability one) the maximal Lyapunov exponent. 
In practice, there will be fluctuations because of many effects, which are discussed in detail in 
[Kantz, 1994| . Based on this understanding, one can derive a robust consistent and unbiased 
estimator for the maximal Lyapunov exponent. One computes 



S{e,m,t) = (in [ —— } \sn+t - Snf+t\ ] ) ■ (19) 





If S{e, m, t) exhibits a linear increase with identical slope for all m larger than some niQ and for 
a reasonable range of e, then this slope can be taken as an estimate of the maximal exponent 
Al. 



13 



The formula is implemented using the Euclidean norm. Apart from parameters characterizing 
the embedding, the initial neighborhood size e is of relevance: The smaller e, the large the linear 
range of S, if there is one. Obviously, noise and the finite number of data points limit e from 
below. It is not always necessary to extend the average in (jl9p over the whole available data, 
reasonable averages can be obtained already with a few hundred reference points s„. If some of 
the reference points have very few neighbors, the corresponding inner sum in ()19p is dominated 
by fluctuations. Therefore one may choose to exclude those reference points which have less 
than, say, ten neighbors. However, discretion has to be applied with this parameter since it may 
introduce a bias against sparsely populated regions. This could in theory affect the estimated 
exponents due to multifractality. Like other quantities, Lyapunov estimates may be affected 
by serial correlations between reference points and neighbors. Therefore, a minimum time for 
|n — n'\ can and should be specified here as well. 




t [section <:rosi3ingsj I [flow samples] 



Figure 5: Estimating the maximal Lyapunov exponent of the C02 laser data (adapted from 
[Hegger et al., 1999| ). Left: results for the Poincare map data, where the average time interval 
Tav is 52.2 samples of the flow, and the straight line indicates A = 0.38. Right: Lyapunov 
exponents determined directly from the flow data. The straight line has slope A = 0.007. 
In good approximation, Xmap = ^fiowTav Here, the time window w to suppress correlated 
neighbors has been set to 1000, and the delay time was 6 units. 

For example, the data underlying the top panel of Figure [5] are the values of the maxima 
of the CO2 laser data. Since this laser exhibits low dimensional chaos with a reasonable noise 
level, we observe a clear linear increase in this semi-logarithmic plot, reflecting the exponential 
divergence of nearby trajectories. The exponent is A ~ 0.38 per iteration (map data!), or, when 
introducing the average time interval, 0.007 per fis [Hegger et al, 1999| . 



2.5.2 The Lyapunov Spectrum 

The computation of the full Lyapunov spectrum requires considerably more effort than just the 
maximal exponent. An essential ingredient is some estimate of the local Jacobians, i.e. of the 
linearized dynamics, which rules the growth of infinitesimal perturbations. One either finds it 
from direct fits of local linear models of the type 

such that the first row of the Jacobian is the vector o„, and 

{J)ij = 5i-ij for i = 2, . . . , m, 
where m is the embedding dimension. The a„, is given by the least squares minimization 

I 
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where {si} is the set of neighbors of s„ [Hegger et al., 1999| . Or one constructs a global nonlinear 
model and computes its local Jacobians by taking derivatives. In both cases, one multiplies 
the Jacobians one by one, following the trajectory, to as many different vectors Uk in tangent 
space as one wants to compute Lyapunov exponents. Every few steps, one applies a Gram- 
Schmidt orthonormalization procedure to the set of Uk , and accumulates the logarithms of their 
rescaling factors. Their average, in the order of the Gram-Schmidt procedure, give the Lyapunov 
exponents in descending order. Apart from the problem of spurious exponents, this method 
contains some other pitfalls: It assumes that there exist well defined Jacobians, and does not 
test for their relevance. In particular, when attractors are thin in the embedding space, some 
(or all) of the local Jacobians might be estimated very badly. Then the whole product can 
suffer from these bad estimates and the exponents are correspondingly wrong. Thus the global 
nonlinear approach can be superior, if a modelling has been successful. 

The computation of the first part of the Lyapunov spectrum allows for some interesting 
cross-checks. It was conjectured jKaplan &: Yorke, 1987| , and is found to be correct in most 
physical situations, that the Lyapunov spectrum and the fractal dimension of an attractor are 
closely related. If the expanding and least contracting directions in space are continuously 
filled and only one partial dimension is fractal, then one can ask for the dimensionality of a 
(fractal) volume such that it is invariant, i.e. such that the sum of the corresponding Lyapunov 
exponents vanishes, where the last one is weighted with the non-integer part of the Kaplan-Yorke 
dimension: 



|Afc+l| 

where k is the maximum integer such that the sum of the k largest exponents is still non-negative. 
Dky is conjectured to coincide with the information dimension. 

The so-called Pesin identity is valid under the same assumptions and allows to compute the 
Kolmogorov- Sinai entropy: 

m 

hKS = ^0(Aj)Aj, 

1=1 

where Q is the Heaviside step function. 

2.6 Dimensions and Entropies of Time Series Data 

Solutions of dissipative dynamical systems cannot fill a volume of the phase space, since dissi- 
pation is synonymous with a contraction of volume elements under the action of the equations 
of motion. Instead, trajectories are confined to lower dimensional subsets which have measure 
zero in the phase space. These subsets can be extremely complicated, and frequently they 
possess a fractal structure, which means that they are in a nontrivial way self-similar. Gen- 
eralized dimensions are one class of quantities to characterize this fractality. The Hausdorff 
dimension is, from the mathematical point of view, the most natural concept to characterize 



fractal sets Eckmann & Ruelle, 1985 , whereas the information dimension takes into account 
the relative visitation frequencies and is therefore more attractive for physical systems. Finally, 
for the characterization of measured data, other similar concepts, like the correlation dimension, 
are more useful. One general remark is highly relevant in order to understand the limitations of 
any numerical approach: dimensions characterize a set or an invariant measure whose support 
is the set, whereas any data set contains only a finite number of points representing the set or 
the measure. By definition, the dimension of a finite set of points is zero. When we determine 
the dimension of an attractor numerically, we extrapolate from finite length scales, where the 
statistics we apply is insensitive to the finiteness of the number of data, to the infinitesimal 
scales, where the concept of dimensions is defined. This extrapolation can fail for many reasons 
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which will be partly discussed below. Dimensions are invariant under smooth transformations 
and thus again computable in time delay embedding spaces. 

Entropies are an information theoretical concept to characterize the amount of information 
needed to predict the next measurement with a certain precision. The most popular one is 
the Kolmogorov-Sinai entropy. We will discuss here only the correlation entropy, which can 
be computed in a much more robust way. The occurrence of entropies in a section on dimen- 
sions has to do with the fact that they can be determined both by the same statistical tool 



Hegger et al., 1999 



2.6.1 Correlation Dimension 

Roughly speaking, the idea behind certain quantifiers of dimensions is that the weight p{e) of 
a typical e-ball covering part of the invariant set scales with its diameter like p{e) ~ e^, where 
the value for D depends also on the precise way one defines the weight. Using the square of the 
probability pi to find a point of the set inside the ball, the dimension is called the correlation 



dimension D2, which is computed most efficiently by the correlation sum [Hegger et al., 1999 



1 ^ 



-^pairs ■ , , . 
^ j=mk<j—w 

where Sj are m-dimensional delay vectors, while 

A^pairs = {N -m + l){N -m-w + l)/2 

is the number of pairs of points covered by the sums, @ is the Heaviside step function and w 
will be discussed below. On sufficiently small length scales and when the embedding dimension 
m exceeds the box-dimension of the attractor [Sauer &: Yorke, 1993] , 

C{m,e) cx 

Since one does not know the box-dimension a priori, one checks for convergence of the estimated 
values of D2 in m. 

Fast implementation of the correlation sum have been proposed by several authors. At small 
length scales, the computation of pairs can be done in 0{N log N) or even 0{N) time rather 
than 0{N'^) without loosing any of the precious pairs. However, for intermediate size data sets 
we also need the correlation sum at intermediate length scales where neighbor searching becomes 
expensive. 

2.6.2 Information Dimension 

Another way of attaching weight to e-balls, which is more natural, is the probability pi itself. 
The resulting scaling exponent is called the information dimension Di. Since the Kaplan- 
Yorke dimension is an approximation of Di, the computation of Di through scaling properties 
is a relevant cross-check for highly deterministic data. Di can be computed from a modi- 
fied correlation sum, where, however, unpleasant systematic errors occur. The fixed mass ap- 
proach [Badii Politi, 1985[ circumvents these problems, so that, including finite sample cor- 



rections [Grassberger, 1988 , a rather robust estimator exists. Instead of counting the number 



of points in a ball one asks here for the diameter e which a ball must have to contain a certain 
number k of points when a time series of length is given. Its scaling with k and A^ yields the 
dimension in the limit of small length scales by 

n / ^ r d\ogk/N 
L)i[m) = lim 



k/N~>o d < \oge{k/N) >' 
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Unlike the correlation sum, finite sample corrections are necessary if k is small. Essentially, the 
log of k has to be replaced by the digamma function ^(k). Given m and r, the routine varies k 
and N such that the largest reasonable range of k/N is covered with moderate computational 
effort. This means that for < k/N < K/N (default: K = 100), all N available points are 
searched for neighbors and k is varied. For K/N < k/N < 1, k = K is kept fixed and is 



decreased Hegger et al., 1999 



2.6.3 Entropy estimates 

The correlation dimension characterizes the e dependence of the correlation sum inside the 
scaling range. It is natural to ask what we can learn form its m-dependence, once m is larger 
than Dq. The number of e-neighbors of a delay vector is an estimate of the local probability 
density, and in fact it is a kind of joint probability: All m-components of the neighbor have to be 
similar to those of the actual vector simultaneously. Thus when increasing m, joint probabilities 
covering larger time spans get involved. The scaling of these joint probabilities is related to the 
correlation entropy /12, such that 

C(m,e) «e^2g-rrxh2^ 

As for the scaling in e, also the dependence on m is valid only asymptotically for large m, which 
one will not reach due to the lack of data points. So one will study h2{m) versus m and try 
to extrapolate to large m. The correlation entropy is a lower bound of the Kolmogorov-Sinai 
entropy, which in turn can be estimated by the sum of the positive Lyapunov exponents. 

An alternate means of obtaining these and the other generalized entropies is by a box counting 
approach. Let pi be the probability to find the system state in box i, then the order q entropy 
is defined by the limit of small box size and large m of 



.1 ~ Q-nihq 



To evaluate YliPi o^^^ ^ mesh of boxes in m >> 1 dimensions, economical use of memory 



is necessary: A simple histogram would take (l/e)™" storage Hegger et al, 1999 
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